N = 16;
ntotal = N*N;
GRD = reshape(1:ntotal,N,N);

Kernel = zeros(ntotal, ntotal);
nWidth = N/4;

Source = [1, nWidth; 1 nWidth];
Target = [2*nWidth+1, 3*nWidth; 2*nWidth+1, 3*nWidth];
SourceInd = reshape( GRD(Source(1,1):Source(1,2), ...
  Source(2,1):Source(2,2)), nWidth * nWidth, 1 );
TargetInd = reshape( GRD(Target(1,1):Target(1,2), ...
  Target(2,1):Target(2,2)), nWidth * nWidth, 1 );

A = randn(nWidth*nWidth,10)*randn(10,nWidth*nWidth) + ...
  1d-10 * randn(nWidth*nWidth, nWidth * nWidth );

R1 = randn(nWidth*nWidth, 20);
Q1 = A*R1;
[U,S,T] = svd(Q1);
U = U(:, find( diag(S) > 1d-6 ));
R2 = randn(nWidth*nWidth, 20);
Q2 = transpose(A)*R2;
[V,S,T] = svd(Q2);
V = V(:, find( diag(S) > 1d-6 ));

Sigma = pinv( transpose(R2)*U ) * transpose(R2)*Q1 * ...
  pinv( transpose(V) * R1 );
A1 = U*Sigma*transpose(V);
norm(A-A1) ./ norm(A)

Kernel(SourceInd, TargetInd) = A;
Kernel = Kernel + transpose(Kernel);

[MatDiag, SampleList, RKTree] = ConstructHMatrix( Kernel, 1 );

nRank = size(RKTree{1}.RKBlockU{2});
nRank = nRank(2);
B = reshape(RKTree{1}.RKBlockU{2}, nWidth*nWidth, nRank ) * ...
  transpose(reshape(RKTree{1}.RKBlockV{2}, nWidth * nWidth, nRank));
norm( A - B )
